﻿ 
using System;
using System.Collections.Generic;
using System.Text;

namespace CO2FixCalibration
{
    static class Statistics
    {
        //static NormalDist maxbiom;
        //static NormalDist caidistr;
        static DataSet[] measuredbiomass;
        static double[] likelihood;
        static double[] rss;
        static  double likelihood_overall;
        public static DataSet[] MeasuredBiomass { get { return measuredbiomass; } }
        public static double[] Likelihood { get { return likelihood; } }
        public static double[] RSS { get { return rss; } }
        
        /*
        public static void Initialize(NormalDist CAIdistr, NormalDist MaxBiom, DataSet[] MeasuredBiomass)
        {
            maxbiom = MaxBiom;
            caidistr = CAIdistr;
            measuredbiomass = MeasuredBiomass;
        }
        */
        public static void Initialize(DataSet[] MeasuredBiomass)
        {
            //maxbiom = MaxBiom;
            //caidistr = CAIdistr;
            measuredbiomass = MeasuredBiomass;

        }
        private static DataSet AggregatedData(DataSet[] AllData)
        {
            DataSet Aggregate = new DataSet("aggr");
            foreach (DataSet d in AllData)
            {
                foreach (int t in d.Time)
                {
                    Aggregate.Add(t, d.GetValueAtTime(t));
                }
            }
            return Aggregate;
        }
        
        public static double Likelihood_biomass(DataSet[] ModeledBiomass, DataSet[] ModeledCAI)
        {
            DataSet SumMeasBiomass = AggregatedData(MeasuredBiomass);
            DataSet SumModelBiomass = AggregatedData(ModeledBiomass);
            double LikelihoodMeas = CalculateLikelihood(SumMeasBiomass, SumModelBiomass);
            return LikelihoodMeas;
        }
        public static double Likelihood_overall(DataSet[] ModeledBiomass, DataSet[] ModeledCAI)
        {
            /*
            rss = new double[MeasuredBiomass.Length];
            likelihood = new double[ModeledBiomass.Length];
            double[] RSS = new double[ReadMeasurements.MeasurementsPerSite.Length];
            
            for (int s = 0; s < ReadMeasurements.MeasurementsPerSite.Length; s++)
            {
                Likelihood[s] = Statistics.CalculateLikelihood(ReadMeasurements.MeasurementsPerSite[s], ModeledBiomass[s]);
                RSS[s] = Statistics.CalculateRSS(ReadMeasurements.MeasurementsPerSite[s], ModeledBiomass[s]);
            }
            */
            DataSet SumMeasBiomass = AggregatedData(MeasuredBiomass);
            DataSet SumModelBiomass = AggregatedData(ModeledBiomass);
            DataSet CAI = AggregatedData(ModeledCAI);

            double LikelihoodMeas = CalculateLikelihood(SumMeasBiomass, SumModelBiomass);
            //double LikelihoodMaxBiom = CalculateLikelihood(SumModelBiomass.MaxValue, maxbiom);
            //double LikelihoodCAI = CalculateLikelihood(CAI.Avg, caidistr);

            //likelihood_overall = LikelihoodMeas * LikelihoodMaxBiom * LikelihoodCAI;

            return LikelihoodMeas;
        }

        private static double CalculateLikelihood(double Value, NormalDist dist)
        { 
            double den = 2* Math.Pow(dist.S2,2);
            double ctr =  Math.Pow((Value - dist.MU), 2);
            double exponent = - (ctr/den);
            double L = 1.0 / (dist.S * Math.Sqrt(2 * Math.PI)) * Math.Exp(exponent);
            return L;
        }
        private static double CalculateLikelihood(DataSet Measured, DataSet Modeled)
        {
            double ctr = CalculateRSS(Measured, Modeled);
            double den = 2* Math.Pow(Measured.Sigma,2);
            double exponent = - (ctr/den);
            return 1.0 / (Measured.Sigma * Math.Sqrt(2 * Math.PI)) * Math.Exp(exponent);
        }
        public  static double CalculateRSS(DataSet Measured, DataSet Modeled)
        {
            int pointsfound = 0;
            return CalculateRSS(Measured, Modeled, out  pointsfound);
        }
        public static double CalculateRSS(DataSet Measured, DataSet Modeled, out int pointsfound)
        {
            pointsfound=0;
            double RSS = 0;
            foreach (int t in Modeled.Time)
            {
                if (Measured.GetValueAtTime(t) > 0)
                {
                    RSS += Math.Pow (Measured.GetValueAtTime(t) - Modeled.GetValueAtTime(t),2);
                    pointsfound++;
                }
            }
            return RSS;
        }
    }
}
